## How do people evaluate foreign aid to “nasty” regimes?
## Tobias Heinrich & Yoshiharu Kobayashi
## British Journal of Political Science
#########################################################


## Run two "basic" regressions

## Load data
############
load("output/Modified MTurk Data.Rdata")
load("output/Modified CCES Data.Rdata")
load("output/Weight vectors.Rdata")


## Estimate Model for non-experimental variables
################################################
all_nonexp_mod <- vector("list", ncol(weight_vec[[1]]))
all_coefs <- c()
for(i in 1:ncol(weight_vec[[1]]))
{
  ## Cluster-bootstrap  
  ## In order to make the comparisons for age nicer, we simply rescale the log(Age) variable.
  w <- weight_vec[[1]][,i]
  mod <- lm(Rating ~ 1 + factor(Income) + factor(Education) + factor(Gender) + factor(RepDem) +
              factor(EconChange) + I(log(Age)-log(18)) + factor(Occupation), 
            data=data, weights=w)
  cluster <- data$ID
  clusters <- unique(cluster)
  nc <- length(clusters)
  Obsno <- split(1:nrow(data), cluster)
  nvars <- length(coef(mod))
  bar <- rep(0, nvars)
  cov <- matrix(0, nvars, nvars)    
  for(iter in 1:n_bs)
  {
    ## Draw clusters
    k <- sample(clusters, nc, replace = TRUE)
    obs <- unlist(Obsno[k])
    
    mod0 <- lm(Rating ~ 1 + factor(Income) + factor(Education) + factor(Gender) + factor(RepDem) +
                 factor(EconChange) + I(log(Age)-log(18)) + factor(Occupation), 
               data=data[obs,], weights=w[obs])
    bar <- bar + coef(mod0)
    cov <- cov + coef(mod0) %*% t(coef(mod0))
  }    
  ## Obtain cluster-adjusted VCOV matrix
  bar <- bar/n_bs
  cov <- (cov - n_bs * bar %*% t(bar))/(n_bs - 1)  
  
  all_nonexp_mod[[i]]$mod <- mod
  all_nonexp_mod[[i]]$VCV <- cov
  
  coef <- data.frame(Name=c("Level 2", "Level 3", "Level 4", "Level 5", 
                            "Level 6", "Level 7", "Level 8", "Level 9", 
                            "Level 10", "Level 11", "Level 12",
                            "Some college", "2 year degree", "4 year degree",
                            "Post graduate", "Female",  "Not very strong Democrat",
                            "Lean Democrat", "Independent", "Lean Republican",
                            "Not very strong Republican", "Strong Republican",
                            "Life got better", "Life stayed the same", "Life got worse", 
                            "Life got much worse", "Age (log)", "Part time employee",
                            "Temp. laid off", "Unemployed", "Retired",  "Permanently disabled",  
                            "Homemaker",  "Student",  "Other employment"),
                     Coefficient=coef(mod)[-1],
                     Lower=coef(mod)[-1] - 1.96 * sqrt(diag(cov))[-1],
                     Upper=coef(mod)[-1] + 1.96 * sqrt(diag(cov))[-1])
  
  coef <- rbind(coef,
                data.frame(Name=c("Level 1", "HS degree or less", "Male",
                                  "Very strong Democract", "Life got much better", "Full time employee"),
                           Coefficient=0, Lower=0, Upper=0))

  
  ## Add age effect
  ind <- which(colnames(cov) == "I(log(Age) - log(18))")
  age_tmp <- matrix(rnorm(3000, coef(mod)["I(log(Age) - log(18))"], sqrt(cov[28, 28])), ncol=3) * (log(matrix(c(18, 50, 70), ncol=3, nrow=1000, byrow=TRUE)) - log(18))
  
  coef <- rbind(coef,
                data.frame(Name=c("Age 18", "Age 50", "Age 70"),
                           Coefficient=colMeans(age_tmp),
                           Lower=colMeans(age_tmp) - 1.96 * colSds(age_tmp),
                           Upper=colMeans(age_tmp) + 1.96 * colSds(age_tmp)))
  coef <- coef[coef$Name != "I(log(Age) - log(18))",]
    
  coef$Header <- 0  
  coef <- rbind(coef, 
                data.frame(Name=c("Income", "Education", "Gender", "Ideology",
                                  "Change in life", "Occupation", "Age"),
                           Coefficient=0, Lower=0, Upper=0, Header=1))

  coef$Name <- as.character(coef$Name)
  coef$Name[coef$Header == 0] <- paste0("   ", coef$Name[coef$Header == 0])

  coef$Formula <- "Complex"
  if(i == 2) coef$Formula <- "Basic"
  if(i == 3) coef$Formula <- "Basic + demographics"
  
  all_coefs <- rbind(all_coefs, coef)
}  
all_coefs$Name <- as.character(all_coefs$Name)
all_coefs <- all_coefs[all_coefs$Name != "Age (log)",]
all_coefs <- all_coefs[all_coefs$Name != "   Age (log)",]

## Saved models
save(all_nonexp_mod, file="output/Saved Nonexperimental Modes.Rdata")


## Prettify variable names
##########################
all_coefs$PrettyName <- factor(all_coefs$Name, 
                    levels=rev(c("Age",
                                 "   Age 18", "   Age 50", "   Age 70", 
                                 "Income", "   Level 1", "   Level 2", "   Level 3", "   Level 4", 
                                 "   Level 5", "   Level 6", "   Level 7", "   Level 8", "   Level 9", 
                                 "   Level 10", "   Level 11", "   Level 12", "Education", 
                                 "   HS degree or less", "   Some college", "   2 year degree", "   4 year degree",
                                 "   Post graduate", "Gender", "   Male", "   Female", 
                                 "Ideology", "   Very strong Democract", "   Not very strong Democrat",
                                 "   Lean Democrat", "   Independent", "   Lean Republican", "   Not very strong Republican",
                                 "   Strong Republican", "Change in life", "   Life got much better", "   Life got better",
                                 "   Life stayed the same", "   Life got worse", "   Life got much worse",
                                 "Age (log)", "Occupation", "   Full time employee", "   Part time employee",
                                 "   Temp. laid off", "   Unemployed", "   Retired", "   Permanently disabled",
                                 "   Homemaker", "   Student", "   Other employment")))


## Graph coefficients for non-experimental results
##################################################
g <- ggplot(data=all_coefs, aes(x=PrettyName, y=Coefficient))
g <- g + geom_point(colour=NA, fill=NA)+ coord_flip() + xlab("") + theme_bw()
g <- g + geom_hline(yintercept=0)
g <- g + facet_grid(~ Formula) + ylab("Effect")
g <- g + theme(strip.text.x=element_text(hjust=0, vjust=.2, size=rel(1.23)),
               strip.text.y=element_text(hjust=0, size=rel(1.23)),
               strip.background=element_blank(),
               axis.text.y = element_text(hjust=0),
               panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                                 linetype="dashed"),
               panel.grid.major.x = element_blank(),
               panel.grid.minor = element_blank(),
               axis.ticks = element_blank(),
               plot.title = element_text(size=rel(1.8), hjust=0))
g <- g + geom_linerange(data=subset(all_coefs, Header == 0),
                        aes(x=PrettyName, ymin=Lower, ymax=Upper, group=Formula),
                        position=position_dodge(width=.8))
g <- g + geom_point(data=subset(all_coefs, Header == 0),
                    aes(x=PrettyName, y=Coefficient, group=Formula),
                    position=position_dodge(width=.8), size=2.4)
ggsave(file="output/figures/NonExpRegression-Long.pdf", width=12, height=8)

